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The  resulting  solution  spectrum  contains  values  which  indicate  the  growth  rates 
of  the  Toilmien-Schlichting  and  flow-induced  surface  instabilities.  It  is  shown  that 
the  Toilmien-Schlichting  instability  is  most  sensitive  to  changes  in  the  surface  prop¬ 
erties.  Previously  it  has  been  suggested  that  an  attempt  to  stabilize  one  class  of 
instability  tends  to  destabilize  the  other  class.  It  is  shown  that  varying  the  surface 
properties  can  reduce  the  growth  rate  of  the  Toilmien-Schlichting  instability  but 
has  little  effect  on  the  flow-induced  surface  instability. 

The  surface  properties  are  '“optimized"  using  a  minimization  algorithm.  It  is 
found  that  appropriate  surface  properties  lead  to  a  decrease  in  the  growth  rates 
of  the  flow  instability^  Although  this  approach  may  be  used  it  is  more  expensive 
computationally  than  a  simple  property  variation  approach. 

The  simple  mechanical  model  for  the  compliant  surface  mayJ>#-represented  by 
an  elastic  plate  over  spring-rigid  supports.  The  functional  relationship  between  the 
flexural  rigidity,  thickness,  and  modulus  of  elasticity  of  the  plate  provides  a  means 
to  vary  the  properties  and  determine  the  effect  on  the  instabilities.  It  is  found  that 
by  keeping  the  flexural  rigidity  essentiaiiy  constant  and  simultaneously  increasing 
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the  plate  thickness  and  decreasing  the  modulus  of  elasticity  a  decrease  in  the  growth 
rate  of  the  Toilmien-Schlichting  instability  is  obtained.  Alternatively,  by  keeping 
the  plate  thickness  and  modulus  of  elasticity  essentially  constant  and  decreasing  the 
flexural  rigidity  a  decrease  in  the  growth  rate  of  the  Toilmien-Schlichting  instability 

results.  Throughout  this  analysis  little  variation  is  found  to  occur  in  the  growth 
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rates  of  the  flow-induced  surface  instability. 

Finally,  the  angle  between  the  rigid  support-arm  and  the  horizontal  in  the 
mechanical  surface  model  is  varied  while  holding  the  surface  properties  constant. 

It  is  shown  that  an  angle  choice  of  between  0  and  50  may  significantly  decrease  the 
growth  rate  of  the  Toilmien-Schlichting  instability. 
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A  spatial  stability  analysis  is  performed  for  the  boundary  layer  over  a  non¬ 
isotropic  compliant  surface.  A  simple  mechanical  model  is  used  for  the  surface. 
Surface  properties  which  may  lead  to  boundary  layer  stabilization  are  determined. 
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of  the  flow  instability.  Although  this  approach  may  be  used  it  is  more  expensive 
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by  keeping  the  flexural  rigidity  essentially  constant  and  simultaneously  increasing 
the  plate  thickness  and  decreasing  the  modulus  of  elasticity  a  decrease  in  the  growth 
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the  plate  thickness  and  modulus  of  elasticity  essentially  constant  and  decreasing  the 
flexural  rigidity  a  decrease  in  the  growth  rate  of  the  Tollmien-Schlichting  instability 
results.  Throughout  this  analysis  little  variation  is  found  to  occur  in  the  growth 
rates  of  the  flow-induced  surface  instability. 

Finally,  the  angle  between  the  rigid  support-arm  and  the  horizontal  in  the 
mechanical  surface  model  is  varied  while  holding  the  surface  properties  constant. 
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CHAPTER  1 


INTRODUCTION 

This  thesis  is  devoted  to  identifying  non-isotropic  compliant  surface  properties 
which  produce  a  delay  in  the  transition  to  turbulence  for  hydrodynamic  applications. 
This  involves  using  a  simple  mechanical  model  for  the  surface.  A  disturbance  is 
introduced  in  the  boundary  layer  in  the  form  of  a  travelling  wave.  The  streamwise 
wavenumber  of  the  disturbance  becomes  the  eigenvalue  for  the  nonlinear  eigenvalue 
problem  formed.  A  measure  of  the  instability  growth  rates  is  found  in  the  solution 
spectrum.  While  the  surface  properties  are  varied  the  least  damped  wavenumber 
is  tracked  to  indicate  the  effect  felt  by  the  instabilities.  It  may  be  shown  that  a 
proper  combination  of  surface  properties  can  lead  to  boundary  layer  stabilization. 

The  transition  of  boundary  layers  from  laminar  to  turbulent  is  due  to  insta¬ 
bilities  that  develop  in  the  boundary  layer.  For  low  Reynolds  number  flows,  the 
viscosity  is  dominant  and  provides  a  means  to  damp-out  the  instability.  As  the 
Reynolds  number  increases  the  natural  damping  becomes  insufficient  to  maintain 
laminar  flow.  Waves  buildup  and  eventually  turn  turbulent.  With  the  onset  of 
turbulence  the  boundary  layer  thickens  and  drag  and  noise  levels  increase.  In  or¬ 
der  to  delay  this  effect  it  is  desirable  to  introduce  a  passive,  “artificial”  damping 
mechanism.  This  may  be  accomplished  by  modifying  the  surface  in  contact  with 
the  boundary  layer. 

A  major  incentive  for  using  a  surface  other  than  a  rigid  wall  was  brought  about 
by  experiments  performed  by  Kramer  jl.2l  in  1960.  By  using  a  rubber  coating  on 
a  rigid  plate,  he  obtained  drag  reductions.  From  his  experiments,  he  concluded: 
(1)  the  surface  induced  artificial  damping  is  a  means  for  boundary  layer  stabi¬ 
lization:  (2)  the  dimensions  and  properties  of  the  elastic  coating  for  an  average 
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Reynolds  number  and  speed  may  be  obtained  through  a  simplified  theory  of  dis¬ 
tributed  damping;  (3)  up  to  a  60  percent  drag  reduction  was  realized  for  the  coated 
surface  compared  to  an  uncoated  identically  shaped  model;  (4)  laminar  recovery  is 
possible  behind  surface  imperfections  which  would  normally  lead  to  transition;  (5) 
no  performance  losses  occurred  due  to  water  impurities;  and  (6)  as  the  Reynolds 
number  increases  the  effect  of  damping  should  increase  due  to  the  improved  contact 
between  the  boundary  layer  and  wetted  surface.  Much  skepticism  has  mounted 
in  reference  to  the  results  of  Kramer  since  experimental  duplication  has  yet  to  be 
realized. 

Important  understandings  of  the  instabilities  occurring  over  a  flexible  surface 
have  been  brought  about  by  the  contributions  of  Benjamin  [3].  His  classification 
of  disturbances  over  a  flexible  surface  was  due  in  part  to  a  stability  discussion  by 
Lin  [4,5,6]  for  two-dimensional  parallel  flows  and  the  analysis  by  Miles  [7,8,9,10]  on 
surface  wave  generation  by  shear  flows.  Landahl  [  1 1  ]  and  Benjamin  [12,13]  further 
identify  distinct  characteristics  which  seperate  the  modes  of  instability  into  three 
classes:  Class  A,  Class  B,  and  Class  C.  The  Class  A  instability  is  realizable  in  the 
presence  of  viscosity  and  is  essentially  a  Tollmien-Schlichting  instability  modified 
by  the  flexible  surface.  The  waves  are  associated  with  a  decrease  of  the  total  kinetic 
energy  of  the  fluid  and  elastic  energy  of  the  wall.  Dissipation  serves  to  increase 
the  wave  amplitude  to  compensate  for  the  energy  loss.  The  waves  are  identified 
as  having  a  speed  less  than  the  velocity  of  the  free  surface  waves  as  was  discussed 
by  Grosch  and  Salwen  (14).  A  Class  B  instability  may  occur  irrespective  of  the 
presence  of  viscosity  and  is  presumed  similiar  to  waves  induced  by  wind  over  water 
surfaces.  Dissipation  in  the  wall  tends  to  stabilize  the  wave.  The  instability  may  be 
recognized  by  a  speed  greater  than  the  free  surface  wave.  And  a  Class  C  instability  is 
realized  where  the  effective  stiffness  of  the  panel  is  too  low  to  withstand  the  pressure 


m 
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forces  induced  on  the  flexible  wall.  This  instability  is  more  commonly  referred  to  as 
a  Kevin-Helmholtz  instability  and  occurrs  when  conservative  hydrodynamic  forces 
cause  a  unidirectional  transfer  of  energy  to  the  solid. 

Grosskreutz  [15,16]  introduced  a  new  approach  in  1971  which  focused  on  the 
control  of  boundary  layer  stabilization  by  the  use  of  non-isotropic  compliant  coat¬ 
ings.  His  experiments  show  that  compliant  coatings  may  lead  to  an  increase  or 
decrease  in  momentum  thickness  which  corresponds  to  an  increase  or  decrease  in 
drag,  respectively.  So  depending  on  the  properties  of  the  compliant  coating  favor¬ 
able  effects  may  be  obtained  or  adverse  effects  may  become  dominant. 

Carpenter  and  Garrad  (17,18)  sought  to  remove  the  skepticism  formed  with 
respect  to  the  isotropic,  Kramer-type  surface  and  expand  on  a  numerical  model 
representing  the  surface.  They  argue  that  a  Kramer  surface  does  have  potential 
for  transition  delay  and  the  reason  skepticism  arose  was  due  to  deficiencies  in  the 
opposing  investigations.  Also,  the  classification  established  by  Benjamin  was  sim¬ 
plified,  or  reclassified,  to  the  following  two  instability  classes  for  a  boundary  layer. 
These  are  the  Tollmien-Schlichting  instability  (TSI)  and  Flow-induced  surface  in¬ 
stabilities  (FISI).  The  FISI  is  basically  the  Class  B  instability  of  Benjamin  and 
Landahi.  They  explain  that  the  Class  C  instability  is  not  found  due  to  boundary 
layer  effects.  In  the  analysis  of  viscous  substrates.  Carpenter  and  Garrad  concluded 
that  a  stabilizing  effect  is  found  for  TSI  in  the  presence  of  a  substrate  and  where 
the  two  modes  coalesce  viscous  substrates  reduce  the  growth  rates  of  instability. 
The  specific  effect  on  boundary  layer  stabilization  bv  a  viscous  substrate  under  a 
Kramer-type  surface  was  investigated  by  Carpenter,  Caster  and  Willis  il9[.  It  was 
found  to  reduce  the  growth  rates  of  the  Tollmien-Schlichting  instability. 

Carpenter  i 20,2 1  i  arrived  at  optimum  surface  properties  for  the  isotropic  case 
which  resulted  in  growth  rates  of  instability  less  than  the  rigid  wall  case.  Carpenter 


[22]  for  the  non-isotropic  case  identified  a  range  of  desirable  surface  properties. 
Carpenter  and  Morris  [23]  for  spatial  wave  growth  and  later  Carpenter  [24]  for 
temporal  wave  growth  observed  growth  rates  of  instability  for  the  non-isotropic 
compliant  surface  less  than  the  rigid  surface.  Morris  [25]  obtained  a  slightly  modified 
model  which  enabled  a  decrease  in  the  nonlinearity  of  the  eigenvalue  problem  of  [23] 
and  [24]  from  an  order  of  six  to  five  in  the  eigenvalue  parameter.  This  model  is 
extended  in  the  present  discussion  to  a  spatial  stability  analysis  to  identify  optimal 
surface  properties  which  may  lead  to  boundary  layer  stabilization. 

The  equations  governing  the  stability  of  flow  over  a  compliant  surface  are  de¬ 
rived  in  Chapter  2.  This  results  in  the  Orr-Sommerfeld  equation  where  the  depen¬ 
dent  variable  is  the  cross-stream  velocity  component  of  an  infinitesimal  disturbance. 
A  simple  mechanical  model  for  the  non-isotropic  compliant  surface  may  be  repre¬ 
sented  by  an  elastic  plate  over  spring-rigid  supports.  The  model  is  chosen  to  mimic 
the  behavior  of  a  compliant  coating  such  as  that  designed  by  Grosskreutz.  The 
coating  consists  of  a  thin  rubber-type  material  covering  stubs  of  a  similiar  material 
and  a  viscous  substrate  fluid  surrounding  the  stubs.  The  equations  governing  the 
motion  of  an  element  of  this  plate  together  with  appropriate  far  field  conditions 
form  the  necessary  boundary  conditions  to  close  the  problem. 

A  spectral  method  approximation  is  introduced  in  Chapter  3  as  a  means  of 
numerical  solution.  The  resulting  matrix  of  equations  forms  a  nonlinear  eigenvalue 
problem  of  degree  five  in  the  eigenvalue  parameter.  Methods  of  solution  are  then 
discussed.  A  model  problem  with  a  known  solution  is  introduced  to  verify  the 
accuracy  of  the  numerical  methods. 

In  Chapter  4.  the  solutions  to  the  eigenvalue  problem  are  discussed  for  the 
model  problem  and  rigid  wall  and  compliant  wall  boundary  layer  cases.  A  compar¬ 
ison  between  the  rigid  and  compliant  case  is  presented  along  with  the  effect  and 


added  cost  arising  due  to  the  addition  of  compliance. 

In  Chapter  5.  the  means  of  obtaining  a  measure  of  the  sensitivity  of  an  eigen¬ 
value  to  surface  property  changes  is  presented  and  an  accuracy  comparison  is  made 
with  a  finite  difference  approximation. 

In  Chapter  6,  the  methods  of  surface  property  optimization  are  formulated. 
The  effect  of  surface  property  selection  for  boundary  layer  stabilization  is  then 
shown.  The  results  are  presented  giving  a  range  of  property  values  which  may  lead 
to  a  delay  of  transition. 
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CHAPTER  2 


P HYSICAL  DESCRIPTION  AND  DERIVATION  OF  PROBLEM 


2.1  Introduction 

Theoretical  investigations  into  the  initial  stages  of  transition  are  based  on  the 
assumption  that  laminar  flows  are  affected  by  small  disturbances.  For  a  boundary 
layer  on  a  solid  body,  these  disturbances  may  physically  be  due  to  wall  roughness 
or  irregularities  in  the  external  flow.  The  question  to  answer  is  whether  the  dis¬ 
turbances  increase  or  decay  in  time  and  space.  If  the  disturbances  decay,  the  main 
flow  is  considered  to  be  stable:  alternatively,  if  the  disturbances  increase  the  flow 
is  considered  to  be  unstable  and  it  is  argued  that  this  then  leads  to  transition  into 
turbulent  flow.  In  this  section  the  theory  of  linear  stability  is  developed  with  the 
object  of  determining  the  flow  conditions  which  may  lead  to  transition. 

2.2  Governing  Equations 

The  problem  to  be  addressed  is  that  of  a  boundary  layer  over  a  smooth,  solid 
surface  immersed  in  an  incompressible,  uniform  flow  with  constant  velocity  and 
pressure.  The  equations  governing  the  flow  are  the  non-linear  Navier-Stokes  equa- 
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where  u.  r,  u>,  and  p  are  instantaneous  now  properties.  In  stability  theory  of  laminar 
Hows  an  infinitesimal  disturbance  is  introduced  on  to  the  laminar  flow  solution. 


Hence,  the  resulting  motion  has  components 


u{x,y,z,t)  =  U(x,y,z,t)  +  u'(x,y,z,t) 
v{x,y,z,t)  =  V(x,y,z,t)  +  v'(x,y,z,t) 
w{x,y,z,t)  -  W(x,y,z,t)  +  w'{x,y,z,t) 
p(x,y,zj)  =  P(x,y,z,t )  +  p'(x,  y,  z,  t), 


(2.3a) 
(2.36) 
(2.3c) 
(2.3  d) 


where  u'.v'.w'  and  p'  are  the  disturbances  and  U,  V,  W  and  P  are  the  laminar  flow 
solutions.  Equations  (2.3)  are  substituted  into  (2.1)  and  (2.2).  It  is  assumed  that 
the  undisturbed  flow  is  a  solution  of  the  Navier-Stokes  equations  and  that  nonlinear 
terms  in  the  disturbance  are  neglected.  The  remaining  terms  result  in  differential 
equations  governing  the  disturbance.  In  boundary  layer  flows  further  stipulations 
may  be  made  which  simplify  the  governing  equations.  The  motion  is  essentially 
two-dimensional  since  Squire  1 26]  showed  that  the  two-dimensional  flow  analysis  is 
more  critical  than  three-dimensional:  the  undisturbed  streamwise  velocity  depends 
on  y  only  (i.e..  i'  =  f  ’(y));  and  the  remaining  two  mean  components,  V'  and  VT.  are 
everywhere  zero.  These  stipulations  describe  a  class  of  flows  known  as  parallel  flows. 
Boundary  layer  flows  may  be  regarded  as  a  good  approximation  to  a  parallel  flow 
because  the  dependence  of  the  velocity,  U .  in  the  streamwise  x-direction  is  much 
smaller  in  comparison  to  the  cross-stream  y-direction.  The  resulting  components  of 
motion  (2.3)  may  be  simplified  to 


u  =  f.'(x.y)  u'(x.yj) 
v  =  v'(x.y.t) 


(2.1a) 

(2.46) 
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w  =  0  (2.4c) 

p  =  P(x)  +  p(x,y,t).  (2  Ad) 


By  substituting  (2.4)  into  (2.1)  and  (2.2),  the  resulting  equations  describe  the 
disturbance  in  a  boundary  layer. 
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dx  + 
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ed  that  far  from  the  wall  in  the  cross-stream  direction 

the  distur- 

bances  vanish. 


u',v'  ,p'  — ►  0  as  y  — *  oo 


(2.6) 


This  assumption  is  necessary  to  satisfy  the  physical  condition  and  is  suitable  for 
securing  boundary  conditions  for  the  resulting  boundary-value  problem  as  will  later 
be  shown. 

The  disturbance  is  assumed  to  be  a  wave  which  propagates  in  the  x-direction. 
The  stream  function  representing  a  single  oscillation  of  the  disturbance  is  assumed 
to  be  of  the  form 


t£’(x,y,f)  =  U^6'<t>  (y)el^z-ut\  (2.7) 

where  the  wave-length  of  the  disturbance  is  A  —  2 n/a  and  the  frequency  of  the  dis¬ 
turbance  is  jj.  The  nondimensional  distribution.  <D,  is  dependent  on  y  only  since  the 


mean  flow  depends  on  y  only.  The  components  of  the  velocity  perturbation  which 
are  obtained  from  (2.7)  may  be  defined  as  partial  derivatives  of  the  streamfunction 
and  given  as 


=  «V(y)e,(oi-wt) 
dy 

(2.8  a) 

v  =  -~  =  -iad>el^z-^\ 
ox 

(2.86) 

where  the  hat  represents  nondimensionalized  disturbances  and  primes  denote  deriva¬ 
tives  with  respect  to  y.  Eliminating  the  pressure  from  (2.5)  and  substituting  (2.8) 
into  (2.5).  a  fourth-order,  ordinary  differential  equation  results  for  the  cross-stream 
velocity  disturbance.  This  is  given  by 

01U  a(y)v"  -  6(y)t-  =  0.  (2.9) 

where  __ 

a(y)  =  -  iR(aU(y)  -  CJ)  -  2a2 

6(y)  =  i  Ra2  (ciU  (y)  -  £)  +  iRaU  (y)  +  a4. 

This  equation  has  nonconstant  coefficients  and  is  commonly  referred  to  as  the  Orr- 

Sommerfeld  equation  which  is  the  stability  equation  for  small  disturbances  in  lam¬ 
inar  flows.  The  equation  has  been  nondimensionalized  with  the  boundary  layer 
displacement  thickness,  b' ,  the  free-stream  velocity,  U0 c,  and  density,  pa.  The 
Reynolds  number  is  given  by 

R  =  — — .  (2.10) 

u 

With  equation  (2.9).  four  appropriate  conditions  are  required  to  obtain  a  so¬ 
lution  for  the  disturbance.  From  (2.6)  where  the  disturbances  vanish  as  infinity  is 
approached  in  the  cross-stream  direction  two  boundary  conditions  result. 


«%).  »*'(y)  —  9 


as  y  —  oc 


(2.i  n 
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In  the  section  following,  the  remaining  two  boundary  conditions  necessary  to 
solve  (2.9)  will  be  obtained.  These  are  the  equations  describing  the  disturbance  at 
the  compliant  surface. 


2.3  Compliant  Boundary  Conditions: 

A  simple  mechanical  model  for  the  non-isotropic  compliant  surface  may  be 
obtained  from  Morris  [25] .  This  is  a  revised  formulation  of  Carpenter  and  Morris 
f 23]  and  Carpenter  [24].  The  concepts  of  the  model  and  derivation  of  the  equations 
of  motion  for  the  disturbance  at  the  surface  follow  and  conclude  with  the  desired 
boundary  conditions.  The  model  consists  of  a  thin,  elastic  plate  supported  by 
hinged  and  sprung  rigid  members  inclined  to  the  horizontal  at  an  angle  0  when  in 
equilibrium.  A  sketch  showing  the  model  is  given  in  Figure  (2.1).  The  motion  of 
the  surface  is  treated  such  that  each  element  of  the  plate  oscillates  in  a  pendulum 
like  motion  at  the  end  of  its  rigid  member.  In  equilibrium,  the  rigid  members  are 
assumed  at  rest.  The  distance  between  each  member  is  assumed  smaller  than  the 
wave-length  of  a  disturbance  normal  to  the  rigid  member.  An  equation  of  motion 
for  a  surface  element  is  desired  which  satisfies  the  constraint  that  the  total  force 
acting  on  the  surface  by  the  mechanical  forces  is  equal  to  the  forces  caused  by  the 
external  fluid  motion  on  the  surface.  Such  an  equation  may  be  given  by 


d2(C60)  d4r]  .  d2C 

Pmb — ^  B— — jco.sO  + 1\  C  60  —  Eb^-^sinO 


dv 


dx 4 


ax‘ 


—  -  r/ cosO  4-  o'cosO  +  r'sinO. 


(2.12) 


where  the  terms  on  the  left  hand  side  of  (2.12)  refer  to  mechanical  forces  and  the 
terms  on  the  right  refer  to  fluid  motion  forces  due  to  viscosity  and  pressure.  For 
the  case  of  an  isotropic  surface,  viscous  interaction  on  the  the  right  side  is  neglected 


and  0  -  ■  0. 
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The  physical  meaning  of  each  term  is  given  respectively  as 
(1.)  rate  of  change  of  momentum  of  the  surface  element 


(2.)  resistance  due  to  bending  stiffness  of  plate 
(3.)  resistance  due  to  spring  stiffness 

(4.)  tension  force  induced  by  relative  motion  of  adjacent  rigid  members 
(5.)  force  due  to  dynamic  pressure  fluctuations 
(6.)  force  due  to  viscous  normal  stress  fluctuations 
(7.)  force  due  to  viscous  shear  stress  fluctuations 
The  variables  in  (2.12)  may  be  defined  as:  i  and  y  are  the  streamwise  and  cross¬ 
stream  coordinates:  <;  and  rj  correspond  to  the  streamwise  and  cross-stream  surface 
displacements:  66  is  the  angular  displacement  of  the  element  relative  to  equilibrium; 
i  is  the  rigid  member  length;  pm  and  b  are  the  plate  density  and  thickness  respec¬ 
tively;  B  and  E  are  the  flexural  rigidity  and  modulus  of  elasticity  of  the  plate;  K  is 
the  spring  stiffness;  and  p' .  o',  and  r'  are  'he  pressure,  viscous  normal  stress,  and 
viscous  shear  stress  fluctuations  on  the  plate  respectively. 

The  necessary  equations  of  motion  for  the  surface  element  are  coupled  by  a 
relationship  between  the  normal  and  tangential  motions  (q.C)  with  the  angular 
displacement.  66.  This  relationship  may  be  given  by 

<;  =  i  66sin6  and  y  =  f  66cos6,  (2-13) 

or 

t  66  —  y / cos6  and  <;  =  ytan6.  (2.14) 

The  normal  displacement  of  the  surface  is  assumed  to  take  the  form 

y  =  6‘7je,(*z~ujn.  (2.15) 


The  continuity  equation  for  the  normal  motion  at  the  surface  implies 


dr\  drl) 
dt  dx 


=  v  at  y  =  rj, 


(2.16) 


or  linearized:  rju  =  tv. 


An  alternative  form  results  by  letting  7  =  t  60  be  represented  in  a  similiar  manner 
as  the  normal  displacement  (2.15).  This  results  in 


v  —  —iujcos0~) . 


(2.17) 


The  continuity  equation  for  the  tangential  motion  of  the  surface  element  implies 


,  d$  d') 1 

U  +  u  =  —  =  — - sind  at  y  =  n, 
dt  dt  y 


(2.18) 


or  linearized:  U'rj  +  u'  —  ~}sinO. 


The  linearization  occurs  at  y  =  0.  Substituting  appropriately  (2.14),  (2.16)  and 


u  —  -v/ia  into  (2.18)  yields 


a(U  (0)cos$  +  iu>sin8)v  +  ujcosdv'  —  0. 


(2.19) 


Equations  (2.12),  (2.16),  and  (2.19)  are  the  resulting  equations  of  motion  and 
boundary  conditions.  Since  the  normal  and  tangential  surface  motion  is  coupled 
one  equation  may  be  eliminated  algebraically. 

Equation  (2.12)  in  nondimensional  form  appears  as 
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The  normal  viscous  stress  perturbation  is  given  by 


or  nondimensionally 


The  viscous  shear  stress  perturbation  is  given  by 


(2.21a) 


(2.21b) 


or  nondimensionally 
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(2.22a) 


(2.226) 


The  pressure  perturbation  is  found  from  the  linearized  cross-stream  component  of 
the  Navier-Stokes  equations 


auv1  +  a2U  v  +  a2Uv'  -  -ia2p  -  —(a2v'  -  v'"),  (2.23a) 

R 

or 
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R  ~  0,2  +  iuRv  +  laRU  0). 


(2.236) 


From  the  continuity  relations  the  following  may  be  defined 


iv1 


1  =  —=T 


(aU  cosO  +  iausind) 


(2.24) 


Substituting  (2.21),  (2.22),  (2.23),  and  (2.24)  into  (2.20)  and  collecting  terms  with 
similiar  powers  of  a,  the  following  boundary  condition  results. 

a5  [Cbcos20i),(O)]  +  a3  [Crsm20t)'(O)] 

_ ^  cosfl 

+a2  \(2uisin6  —  3iU  (O)cos0) — —  v^O)] 

R 


4 -a[(Cff  —  u!2Cm)v'(0 )  +  {U  (0)cos<?  4-  iusind)  — —  {/'(O)] 
+i(U  (0 )cosd  +  iusind) — -—v"'(0)  -  iu2 sinOcosOv' (0)  =  0 

it 


The  final  boundary  condition  is  given  by  equation  (2.19)  which  is 


(2.25) 


a(U  (O)cos0  4-  ujsm0)t)(O)  4-  wcos0i/(O)  =  0. 


The  equations  of  motion  governing  the  stability  of  the  flow  over  a  non-isotropic 
compliant  surface  for  the  cross-stream  velocity  component  of  the  disturbance  have 
been  derived.  The  means  of  obtaining  the  solution  will  be  described  in  the  next 
chapter.  For  convenience,  the  overbar  on  U  representing  the  Blasius  solution  will 
hereafter  be  neglected. 
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CHAPTER  3 


NONLINEAR  EIGENVALUE  PROBLEM 


3.1  Introduction 

In  the  previous  chapter  a  detailed  analysis  was  performed  to  arrive  at  a  math¬ 
ematical  representation  of  the  physical  problem,  namely  the  stability  analysis  of 
flow  over  a  compliant  surface.  The  present  chapter  will  formulate  the  equations 
into  a  form  suitable  for  obtaining  a  numerical  solution  and  describe  the  methods  of 
solution. 

3.2  The  Orr-Sommerfeld  Problem 

In  hydrodynamic  stability  theory,  the  Orr-Sommerfeld  equation  (2.9)  governs 
the  normal  velocity  component  of  the  disturbance  imposed  on  the  flow.  The  so¬ 
lution  will  give  the  characteristic  instabilities  present  for  a  time  or  space  varying 
analysis.  The  problem  at  hand  varies  in  space  and  is  referred  to  as  a  spatial  stability 
problem.  The  wavenumber,  a  ,  is  complex  and  is  taken  to  be  the  unknown  eigen¬ 
value.  A  negative  imaginary  part  of  a  indicates  that  the  solution  is  growing  in  the 
streamwise  direction.  This  is  an  indication  of  an  instability  growth  present  in  the 
flow.  A  zero  imaginary  part  suggests  that  the  solution  is  neutral.  And  a  positive 
imaginary  part  suggests  that  the  solution  decays  in  the  streamwise  direction.  The 
frequency,  u7,  and  Reynolds  number,  R,  are  both  real  and  specified.  The  alternative 
problem  in  stability  theory  which  will  not  be  solved  in  this  study  is  time  varying, 
or  temporal.  In  such,  a  and  R  are  both  real  and  specified  while  the  frequency 
is  complex  and  becomes  the  unknown  eigenvalue  of  interest.  The  frequency  and 
wavenumber  are  related  and  together  form  the  phase  velocity,  c.  which  is  defined 


The  problem  of  boundary  layer  flow  over  a  flat  plate  is  an  eigenvalue  problem  in 
a  and  is  said  to  be  nonlinear  to  a  degree  of  four  in  the  eigenparameter.  The  problem 
at  hand  is  nonlinear  to  a  degree  of  five  where  the  added  degree  of  nonlinearity  is 
introduced  in  the  boundary  condition  (2.25).  The  nonconstant  coefficients  result 
due  to  the  streamwise  component  of  the  Blasius  velocity  profile  the  solution  of 
which  is  only  known  numerically.  So  the  eigenvalue  problem  must  also  be  solved 
numerically.  The  technique  used  to  formulate  the  numerical  approximation  was 
previously  used  by  Bridges  and  Morris  [ 27,28]  and  Bridges  (29)  in  the  solution  of 
the  fourth-order  nonlinear  eigenvalue  problem  for  the  flow  over  a  rigid  surface.  The 
reason  for  such  an  approach  will  be  made  evident. 

The  domain  of  the  equation  is  from  zero  at  the  surface  to  infinity  in  the  cross¬ 
stream  direction.  In  order  to  solve  the  problem  numerically,  the  domain  may  either 
be  truncated  or  transformed  to  some  finite  domain.  Grosch  and  Orszag  [30i  have 
performed  a  study  of  this  subject  and  suggest  an  algebraic  transformation 


2  =  - -  and  y  =  L  - - .  (3. l«.b) 

y  4 -  L  1-2 

where  26  [-1,  4-  ll  and  y  lO.oc).  In  this  analysis  a  value  of  L  2  will  be 
used  and  is  suggested  as  optimum  by  Bridges  for  the  rigid  surface  problem.  The 
corresponding  metric  arrived  at  is 


m(z)  = 


dz 

dy 


A  far  field  condition  (2.7)  in  the  transformed  domain  appears  as 


m 


«*'( 


(1  -I2  </,-•(  2) 

2  L  '  dz 


(3.2) 


13.3) 


•  0  r  J.s  2  •  1 


This  introduces  an  ambiguity  as  to  the  value  of  t)' ( 1 )  since  the  metric  approaches 


zero  as  z  — *  1.  In  order  to  temporarily  avoid  such  a  problem,  a  nondimensional 
dummy  variable,  £(z)  ,  is  introduced  and  defined  as 

c  —  rnv  .  (3.4) 

Making  the  appropriate  substitutions  of  (3.2)  and  (3.4)  into  (2.9),  the  Orr- 
Sommerfeld  equation  in  the  transformed  domain  may  be  written  as 

m(m(m£')')'  -  ma(z)£'  b(z)i v  =  0, 

where 

a(z)  —  -  iR(al'(z)  -  !•)  -  2 a2 
b(z)  =  t  Ra2  (aU  (z)  -  CJ)  +  iaRm(mU'{z))'  +  a4 

with  far  field  conditions 

r(l)  =  s(l)  =  0  (3.6a, b) 

and  compliant  boundary  conditions 


(3.5a) 

(3.5  b) 
(3.5c) 
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0  (3.6c) 


a{cos$U'(-l)  +  iuisin8)v(-l)  +  iocos6£{—  1)  =0,  (3.6  d) 

where  the  nondimensional,  constant  coefficients  were  previously  defined  as 


n  Pmb 

CM  =  - r  > 

Po<5“ 


CB  = 
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PoU>6‘ 3’ 


C  K 


Ki 5* 

P*UZ, 


and  Ct  = 


Eb 


PoU'8- 


The  primes  in  (3.4)-(3.6)  denote  derivatives  with  respect  to  the  transformed  vari¬ 
able.  z,  and  i'[z)  is  the  Blasius  streamwise  velocity  profile  in  the  transformed  do¬ 
main.  A  plot  of  the  Blasius  profile  versus  the  transformed  variable,  z,  is  shown  in 
Figure  (3.1).  In  the  limiting  case  where  C\i  — -  oo,  the  compliant  problem  becomes 
a  rigid  wall  problem. 

A  spectral  approach  known  as  a  finite  Chebvshev  series  expansion  is  sought  for 
the  solution  of  (3.4)-(3.6).  A  spectral  expansion  is  an  approximation  of  an  unknown 
function  by  a  series  of  known  functions  which  satisfy  the  boundary  conditions. 
Gottlieb  and  Orszag  [3 1 1  and  Fox  and  Parker  [32]  discuss  in  detail  the  advantages 
of  such  an  approximation  and  give  various  examples.  Gottlieb  and  Orszag  state 
that  a  Chebvshev  polynomial  expansion  gives  a  good  representation  of  functions 
that  undergo  rapid  changes  in  narrow  boundary  layers.  One  reason  is  that  the 
polynomials  can  resolve  changes  over  distances  of  order  n  ~2  where  n  is  the  number  of 
Chebvshev  series  terms  retained.  Also,  for  the  Chebyshev  series  expansion  the  error 
converges  exponentially  in  comparison  to  finite  difference  methods  which  converge 
algebraically. 

With  this  in  mind,  the  Chebyshev  series  expansion  for  the  disturbance  in  equa¬ 


tions  (3.4)-(3.6)  may  be  given  by 


(3.7) 


°i2)  =  'vnTn{z) 

ri5=0 

N 

£(*)  =  ^T'Cn^)*).  (3-8) 

n^O 

The  prime  on  the  summation  signifies  that  the  leading  term  of  the  series  is  to  be 
halved.  The  Blasius  velocity  profile  is  expanded  in  a  similiar  series. 


OO 

U{z)  =  Y,'UnTn(z)  (3.9) 

n=0 

Details  on  how  a  known  function  may  be  represented  by  a  Chebyshev  series  expan¬ 
sion  may  be  found  in  Appendix  A  and  in  Appendix  B  specifically  for  the  Blasius 
solution. 

Due  to  the  properties  of  Chebyshev  polynomials  it  is  convenient  to  pose  equa¬ 
tions  (3.4)-(3.6)  in  integral  form.  As  such,  the  following  equations  result. 


mv  — 


i  +  e„ 


(3.10) 


and 


P,(+  J  P,(  +  f  f  +  f  f  J  P,(  +  f  f  J  J  f  Wi¬ 


lli 


bv  +  e\  +  e2Z  +  e3~  =  0 


(3.11] 


where  a(z)  and  b(z)  are  given  by  (3.5b)  and  (3.5c),  and 
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P i  (2)  =  -  6 m2m' 
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P2(2)  =  Tm(m')  -4m 
Pz{z)  =  -  (m(mm')')' 


(3.126) 

(3.12c) 

(3.12(f) 


In  connection  with  a  Chebyshev  series  expansion,  the  Tau  method  which  was 
introduced  by  Lanczos  in  1938  will  be  used  to  remove  the  constants  of  integration. 
A  detailed  explanation  may  be  found  in  [31]  and  [33],  but  at  present  only  a  brief 
outline  of  the  method  will  be  presented  without  proofs. 

The  series  approximation  of  a  function,  £(2)  or  v(z),  previously  introduced  in 

(3.8)  and  (3.7)  has  k  additional  terms  added  to  it  where  k  represents  the  number  of 
independent  boundary  constraints  that  must  be  applied  (i.e.,  one  and  three  in  this 
particular  problem).  The  resulting  approximation  is  the  exact  solution  to  a  slightly 
modified  problem.  This  results  in  (V  +  3  unknowns  for  N  +  2  equations  and  one 
boundary  constraint  for  equation  (3.10)  and  (V  4-7  unknown  coefficients  for  N  +  4 
equations  and  three  boundary  constraints  for  equation  (3.11).  Respectively,  the 
equations  involving  the  coefficients  of  T0(z)  and  T0(z),  T\  (2)  and  T2(*)  for  (3.10) 
and  (3.11)  serve  to  determine  the  constants  of  integration  only  and  so  may  be 
disregarded  for  the  present  analysis.  The  added  “tau"  terms  need  not  be  explicitly 
calculated  either.  The  remaining  system  is  composed  of  N  equations  with  one 
boundary  constraint  and  N  —  2  equations  with  three  boundary  constraints.  The 
two  systems  of  equations  when  combined  result  in  a  square  N+ 1  matrix  of  equations 
as  will  be  shown. 

Using  the  Chebyshev  product  and  integral  formulae,  the  series  expansions  (3.7)- 

(3.9)  and  the  metric  (3.2),  represented  by  the  following  Chebyshev  series 


are  substituted  into  equations  (3.10)  and  (3.11).  This  results  in  a  set  of  equations 
with  the  vectors  of  unknown  Chebyshev  coefficients.  {£}  and  {v}  .  Using  the  far 
field  condition.  £(1)  =  0  .  with  equation  (3.10).  the  following  relation  is  found. 

{£}  =  [T\{v}  (3.14) 

{£}  and  {t>}  are  column  vectors  containing  unknown  Chebyshev  coefficients  and  [Tj 
is  a  square  iV  4-  1  matrix.  The  remaining  three  boundary  conditions  with  equation 
(3.11)  give 

-  5  -j  5 

Y^Ck^~k  {v}^-  {(}  =  { 0},  (3.15) 

'-k=0  Lfc  =  0  J 

k*l 

where  [Ck\  and  \Dk\  are  complex  square  matrices  of  order  N+l  which  are  functions 
of  uJ.  R  and  the  compliant  boundary  condition  properties.  The  dummy  vector.  {£}, 
may  be  eliminated  from  (3.15)  by  the  substitution  of  (3.14).  The  following  nonlinear 
eigenvalue  problem  results, 

[05(a)]  {f}  =  {0},  (3.16) 

where 

Dn{a)  =  C0ah  +  C,a4  +  C2cT  +  C352  +  C4a  +  C5.  (3.17) 

This  forms  the  Chebyshev  discretization  of  the  Orr-Sommerfeld  equation  over  a 
compliant  surface. 


3.3  Solution  of  the  Nonlinear  Matrix  Eigenvalue  Problem 

The  eigenvalue  problem  considered  is  nonlinear  in  a  to  the  degree  of  five  where 
the  highest  degree  of  nonlinearity  is  introduced  in  the  boundary  condition.  The 
system.  Ds(a).  may  be  referred  to  as  a  lambda  matrix.  Since  one  boundary  condi¬ 
tion  is  independent  of  a  it  may  be  eliminated  using  appropriate  column  operations; 
thus,  the  problem  is  reduced  to  N  equations  and  N  unknowns,  or  a  system  of  com¬ 
plex  square  matrices  of  order  N.  For  the  solution  of  the  lambda  matrix,  three  global 
methods  and  one  local  refinement  method  will  be  considered.  A  global  method 
is  global  only  in  the  sense  that  an  initial  guess  for  the  eigenvalue  determination  is 
unnecessary.  In  a  local  method  an  initial  guess  is  required.  The  global  methods  are: 
(1)  linearization  by  a  companion  matrix,  (2)  factorization  with  Bernoulli  iteration 
to  obtain  a  subset  of  the  spectrum,  and  (3)  factorization  with  Traub  iteration  to 
obtain  a  subset  of  the  spectrum.  The  local  method  is  a  refinement  of  Newton’s 
method  derived  by  Lancaster  1 34 1  for  a  single  eigenvalue. 

The  companion  matrix  method  has  been  used  for  the  Orr-Sommerfeld  problem 
by  Benney  and  Orszag  135!.  Bridges  and  Morris  1331  and  Gohberg,  Lancaster  and 
Rodman  (36l  discuss  both  the  companion  matrix  method  and  factorization.  From 
such,  the  analysis  is  extended  to  the  larger  system  at  hand.  The  companion  matrix 
is  a  linearization  of  the  lambda  matrix  and  therefore  is  of  a  larger  order.  If  m  is 
the  order  of  the  matrix  system.  D$(a),  then  the  order  of  the  companion  matrix  is 
5m.  When  a  differential  equation  is  formulated  as  a  matrix  problem,  it  takes  on 
the  form  of 


Ax  =  A  Bx 


(3.19) 


where  A  is  the  eigenvalue  and  x  represents  the  eigenvector.  Eigenvalue  determina- 


tion  is  found  by  the  condition 


Det\A-XB\=0  (3.20) 

A  similiar  construction  for  the  present  problem  yields 


r-Ci 

-c2 

-C3 

-Ca 

— c5 ' 

rC0 

0 

0 

0 

O' 

'  a4a 

i 

0 

0 

0 

0 

0 

I 

0 

0 

0 

a3  a 

0 

/ 

0 

0 

0 

-  a 

0 

0 

I 

0 

0 

a2  a 

0 

0 

I 

0 

0 

0 

0 

0 

/ 

0 

aa 

.  0 

0 

0 

I 

0  J 

.  0 

0 

0 

0 

/  J 

a  , 

(3.21) 

Referring  to  (3.20),  if  B  is  invertible  a  more  efficient  and  equivalent  form  is 

Det\B~xA  -  A/|  =  0.  (3.22) 


The  leading  coefficient  matrix,  [C0],  is  singular  since  the  only  entries  are  introduced 
in  the  compliant  boundary  condition  as 
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(3.23) 


To  remove  the  singularity  in  [C0j  an  algebraic  transformation  is  introduced 


A  = 


1 


a  -  s 


(3.24) 


where  s  is  a  real  constant  taken  in  this  analysis  to  be  (u/0.35).  The  problem  may 
now  be  cast  in  the  form  of  (3.22)  giving 
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(3.25) 


The  eigenvalues  of  (3.25)  may  be  obtained  using  the  efficient  QR  algorithm. 

The  second  method  is  derived  by  a  matrix  equivalent  to  synthetic  division  to 
compute  the  dominant  solvent.  After  applying  the  algebraic  transformation  (3.25), 
the  following  results 


D5(X)  =  {Q4(A)}(A7  -  Y)  +  Rr.  (3.26) 

where 

<?4(A)  =cx  +  (C0Y  +  C  i)A3  +  {C0Y2  +  c,y  +  c2)a2 
+(c0y3  +  c,K2  4-  c2y  +  c3)a 
+(c0y4  +  c,r3  +  c2y2  +  c3y  +  c4) 

and  is  considered  to  be  the  right  quotient  and  /2r  is  the  right  remainder  of  the 
division  of  D 5(A)  by  (XI  -  V').  For  (A/  -  V')  to  be  a  factor  of  Ds(A),  the  remainder. 
Rr,  must  be  set  to  zero.  This  is  given  by 

Rr  =  c0y 5  +  c  1  y 4  +  c2y3  +  c3y2  +  c4y  +  c5  =  o.  (3.27) 

The  square  matrix,  Y,  is  referred  to  as  the  right  solvent.  The  Bernoulli  iteration 
method  will  be  incorporated  to  solve  the  matrix  polynomial  (3.27).  For  such  we 
seek  the  dominant  solvent  which  may  be  obtained  from  the  iteration  formula 


yt+1  =  -C~ 1  (Ci  +  (C2  4-  (C3  -  (C4  4-  C,Y^)Y-l2)Y~\)Yt-1). 


(3.28) 


where  Y0  =  Yx  =  Y2  =  F3  =  0  and  F4  =  -Cn_1C,. 

Upon  convergence,  the  eigenvalues  are  obtained  by  using  the  QR  algorithm. 


The  final  global  method  to  be  considered  was  developed  by  Dennis,  Traub  and 
Weber  [37]  to  compute  a  dominant  solvent.  The  algorithm  is  a  generalization  of 
an  algorithm  for  scalar  polynomials  by  Traub  [38].  The  method  was  discussed  by 
Morris  [25]  for  the  compliant  problem  approaching  the  limiting  case  of  the  rigid 
wall  problem.  The  method  consists  of  two  iterative  steps.  The  first  consists  of 
constructing  the  equivalent  of  the  G-polynomials. 


G0(F)  =  I 

Gn+l(Y)  =  Gn(Y)Y  -r\n)D5(Y), 


(3.29a) 

(3.296) 


where 


Gn(Y)  =  r(,n)r5  +  rin)y4  +  r£n)r3  +  rin)r2  +  r(5n)r  +  r^n).  (3.29c) 


The  second  stage  is  given  by 


Yo  =  (r(/'))(r(Ix'-1)r1  (3.30a) 

and 

Yi+i  =GL(Y,)Glll{Yl)  (3.30b) 

where  L  is  the  final  G-poiynomial  built-up.  The  first  stage  of  the  algorithm  is 
equivalent  to  the  Bernoulli  iteration.  The  second  stage  is  only  linearly  convergent, 
but  the  asymptotic  error  constant  may  be  made  as  small  as  desired  by  increasing 
the  number  of  iterations  of  the  first  stage.  A  subset  of  the  eigenvalue  spectrum  may 
be  obtained  by  using  the  QR  algorithm. 


y.c.V.  /.v; 
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The  final  method  to  be  considered  is  a  locally  convergent  algorithm  which 
requires  a  sufficiently  good  initial  guess  for  the  refinement  of  a  single  eigenvalue. 
The  local  scheme  is  a  refinement  of  Newton's  method  and  has  quadratic  convergence. 
The  method  is  attributed  to  Lancaster  [34 j  and  an  example  of  its  implementation 
may  be  found  in  Bridges  and  Morris  [33].  The  iterative  formula  is  given  by 

5i+i  =5<  -2/(at)/{[/(a,)]2  - /(1)(a,)},  for  i  =  0,1,2, ...  (3.31a) 

where 

/(5.)  =rr{/r1(atp(1)(at)}  (3.316) 

and 

/(l) (<*,•)  =  Tr{D-l(ai)D^(at)  -  (£>_l(a1)D(1)(a,)]2}.  (3.31c) 

Tr{s}  denotes  the  trace  of  matrix  (sj,  D~l  is  the  inverse  of  D  and  D ^  denotes 
the  ith  derivative  of  D  with  respect  to  a.  It  should  be  noted  that  only  one  matrix 
inverse  is  required.  Also,  as  will  be  discussed  in  a  later  chapter,  the  eigenvectors 
necessary  for  the  surface  property  optimization  may  be  conveniently  computed  as 
an  offshoot  of  this  method  making  use  of  the  matrix  operations  already  performed. 

This  concludes  the  outline  of  methods  considered  for  eigenvalue  determination. 
Actual  global  method  comparisons  for  accuracy  and  efficiency  were  not  in  the  main 
context  of  this  investigation.  The  global  schemes  are  necessary  to  determine  a  good 
initial  guess  for  the  least  damped  eigenvalues  for  TSI  and  FISI  to  be  refined  in 
the  more  efficient  local  method.  The  sensitivity  of  the  eigenmodes  to  changes  in 
the  surface  properties  may  then  be  performed,  followed  by  the  optimization  of  the 
surface  properties.  With  respect  to  comparisons  and  applications,  reference  may  be 
made  to  Bridges  and  Morris  [33]  and  Morris  [25]. 
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3.4  A  Model  Eigenvalue  Problem 

When  investigating  a  complex  problem  requiring  numerical  techniques  as  a 
means  to  obtaining  a  solution,  it  is  advantageous  to  devise  a  model  problem  with  an 
exact  solution  which  captures  as  many  of  the  characteristics  of  the  physical  problem 
of  interest  as  possible.  One  chooses  a  model  problem  such  that  the  numerical 
solution  may  be  compared  with  the  known  exact  solution.  The  compliant  surface 
problem  has  many  identifiable  characteristics  most  of  which  may  be  incorporated 
in  the  model. 

The  model  boundary  value  problem  is  given  by 

t<t>"  -  2 au<p'  +  a2 d>  =  0  (3.32a) 

with  boundary  conditions 

o(l)  -  0  (3.326) 

a3<6(  -  1)  +  o'(  -  1)  =  0.  (3.32c) 

where  the  primes  represent  derivatives  with  respect  to  z  [  —  1,4-1],  The  eigenvalue, 
a.  enters  the  boundary  condition  at  a  higher  power  than  in  the  differential  equation 
which  is  similiar  to  the  physical  problem.  A  stiffness  parameter,  e,  which  may  be 
thought  of  as  R~l/2,  multiplies  the  highest  derivative  so  as  to  simulate  the  viscous 
terms  in  the  physical  problem. 

The  exact  solution  of  (3.32)  is  given  by 

<*>(z)  =  (e^z  -  (a3  +  i2)e-^z)/(e^z  +  (a3  +  73)^"*)  (3.33) 

with 

(a3  *  72)e"n  ~  («''’  +  '73 ) e —  ^ 1  =  0  (3.34) 


A  *, 


For  a  numerical  solution  the  equation  may  be  put  in  integral  form,  or 


(3.36) 

<t>(z)  =  ]T'<£nrn(2).  (3.37) 

71  =  0 

By  substituting  (3.37)  into  (3.36)  and  incorporating  the  Tau  method,  the  problem 
may  be  cast  into  a  lambda  matrix  which  is  of  order  three  in  the  eigenvalue  and 
given  by 

Dz  (a)  =  C (jft3  +  C[Qr2  4-  C  2  ot  +  C3 ,  (3.38) 

where  [Cx]  are  complex  square  matrices  of  order  N  +  l.  The  leading  coefficient  ma¬ 
trix,  [C0 1,  is  singular  so  transformation  (3.24)  is  applied.  The  methods  discussed  in 
the  previous  section  apply  in  a  similiar  manner  to  that  of  the  physical  problem. 

The  stage  has  now  been  set  for  solving  the  problem  at  hand  and  the  methods  of 
solution  for  the  nonlinear  eigenvalue  problems  have  been  described.  The  accuracy 
of  the  eigenvalues  and  eigenfunction  for  a  given  number  of  Chebyshev  polynomials 
will  be  tested  for  the  compliant  wall,  the  rigid  wall,  and  the  model  problem  in  the 
next  chapter. 


e<t>  —  2aui 


h+*Jf 


<j>  +  e„  4-  ei2  =  0. 


The  function  is  approximated  by  a  finite  Chebyshev  series 


CHAPTER  4 


NUMERICAL  RESULTS  OF  EIGENVALUE  PROBLEMS 
4.1  Model  Problem 

As  was  mentioned  in  the  previous  chapter,  a  model  boundary-value  problem 
with  characteristics  similiar  to  the  physical  problem  is  used  to  test  the  numerical 
methods.  The  global  methods  of  solution  were  discussed  in  the  previous  chapter 
and  results  from  each  method  are  given  in  Table  (4.1)  for  e  =  1.0  and  u 7  —  0.25. 
As  is  shown,  for  a  small  number  of  Chebyshev  polynomials  the  methods  give  a 
comparatively  similiar  accuracy  for  the  given  number  of  iterations.  The  Bernoulli 
and  Traub  iteration  methods  result  in  only  a  subset  of  the  eigenvalue  spectrum.  As 
can  be  seen  the  third  eigenvalue  is  undetected  by  these  methods  for  N=5.  Since 
only  a  sufficiently  good  initial  guess  is  required  for  an  eigenvalue,  little  else  will 
be  needed  in  the  form  of  demonstration  and  comparison  with  respect  to  the  global 
methods.  A  more  indepth  comparison  of  the  these  methods  and  the  local  method 
may  be  obtained  from  Morris  ! 25] ,  Bridges  and  Morris  1 33] .  Benney  and  Orszag  [ 3 5 1 . 
and  Dennis.  Traub  and  Weber  [37!. 

The  corresponding  eigenfunction  is  obtained  for  the  smallest  eigenvalue  in  Table 
(4.1).  It  is  sufficient  at  present  to  view  the  accuracy  of  the  method  for  a  given 
number  of  Chebyshev  polynomials  (N)  and  an  imposed  stiffness  (c).  In  Figure  (4.1) 
a  plot  for  N  =  5,7,  and  10  with  c  =  1.0.  shows  the  eigenfunction  to  be  somewhat 
independent  of  the  number  of  polynomials.  The  numerical  solution  is  essentially 
indistinguishable  from  the  exact  solution.  The  problem  is  made  stiff  by  requiring 
the  parameter  (  to  be  small.  The  corresponding  eigenfunctions  for  N  -10  and 
(  ~  1.1/ v  100.  and  1  \  500  are  shown  in  Figure  (4.2).  The  numerical  and  exact. 


solutions  again  prove  to  be  indistinguishable.  This  provides  much  encouragement 


=  0 


0.54543587 

0.98403404 

1.41189790 

1.76754890 


N=7 


0.54541156 

0.98403205 

1.41193160 

1.76666540 


N  =  10 


0.54541160 

0.98403009 

1.41191660 

1.76673320 


Exact 


0.54541160 

0.98403009 

1.41191660 

1.76673320 


Bernoulli  iteration(15): 


0.54543587 

0.98408419 

1 .76755 180 


0.54541156 

0.98402771 

1.41134490 

1.76666740 


0.54541160 

0.98407944 

1.41176970 

1.76673320 


0.54541160 

0.98403009 

1.41191660 

1.76673320 


Traub  iteration(5:4): 

N  =  10 


Exact 


0.54543587 

0.54541156 

0.54541160 

0.54541160 

0.98399391 

0.98403206 

0.98402854 

0.98403009 

1.41192940 

1.4119248 

1.41191660 

1.76754910 

1.76666540 

1.76671350  1 

1.76673320 
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for  the  use  of  this  approach  for  the  physical  problem  which  is  stiff. 

4.2  Rigid  Wall  Case 

The  flow  properties  used  are:  the  freestream  velocity  and  density  are  20  m/s 
and  1000  kg/m3,  respectively ;the  Poisson  ratio  is  0.5;and  the  viscosity  is  .001002 
kg/ms.  These  are  obvious  choices  for  the  density  and  viscosity  since  water  can 
be  found  to  have  a  density  and  viscosity  of  1000-1020  kg/m3  and  .001002  kg/ms, 
respectively.  The  compliant  surface  model  is  taken  to  be  a  flexible  plate.  As  the 
mass  of  the  plate  is  increased,  the  characteristics  of  the  problem  become  more 
similiar  to  that  of  a  rigid,  flat  surface.  The  rigid  surface  is  achieved  in  the  limit  as 
the  mass  coefficient,  Cm,  approaches  infinity.  The  solutions  obtained  in  this  limit 
should  coalesce  with  published  results  for  the  solution  of  the  Blasius  velocity  profile 
over  a  flat  plate.  A  common  reference  for  comparison  is  the  neutral  stability  curve. 
For  this  comparison  it  is  sufficient  to  use  only  fifteen  Chebyshev  polynomials  to 
obtain  an  adequate  accuracy.  The  results  are  listed  in  Table  (4.2)  and  a  comparison 
is  made  in  Figure  (4.3)  with  values  from  Jordinson  [39]  and  Van  Stijn  and  Van  De 
Vooren  [40).  As  is  expected  the  results  fall  on  a  common  curve. 

As  with  the  eigenvalue,  the  normalized  eigenfunction  comparison  must  be  made 
with  other  numerical  results  since  the  exact  solution  is  not  known.  Jordinson  [39] 
referred  to  the  case  where  R=998,  u  =  .1122,  and  a  =  (.3086,  -.0057).  Using  N  =  15, 
a  comparison  results  in  nearly  an  exact  fit  as  shown  in  Figure  (4.4).  The  results 
begin  to  deviate  slightly  as  the  distance  from  the  wall  increases.  An  observable 
jump,  or  step,  occurrs  in  Jordinson's  analysis  which  is  not  found  in  the  present 
calculations.  Since  the  function  is  well  behaved  no  physical  explanation  justifies 
such  a  step.  And  finally  a  last  comparison  is  made  to  determine  the  eigenfunction 
sensitivity  to  the  number  of  Chebyshev  polynomials.  Figure  (4.5)  is  an  eigenfunction 


CO 


Table  4.2:  Values  of  R,ar,  and  u>  for  the  neutral  curve  in  the  limit  as  C\i 
which  as  the  compliant  surface  becoming  a  rigid  plate.  (  N  =  15  ). 


R 

olt 

U / 

2200.0 

.3095 

.1010 

1400.0 

.3356 

.1185 

0800.0 

.3557 

.1368 

0520.0 

.3014 

.1193 

0536.5 

.2753 

.1067 

0604.0 

.2406 

.0893 

1364.0 

.1450 

.0433 

-.1 
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comparison  for  N=  10,15,  and  20  with  the  results  by  Jordinson  corresponding  to 
R=336,  u7  =  .1297,  and  a  =  (.3084,  .0079).  As  with  the  model,  the  eigenfunctions 
are  relatively  independent  of  the  number  of  Chebyshev  polynomials  for  the  rigid 
wall  case  to  achieve  sufficiently  accurate  solutions. 


4.3  Compliant  Wall  Case 

The  discussion  in  this  section  will  primarily  be  devoted  to  determining  an 
adequate  number  of  polynomials  required  for  a  sufficiently  accurate,  or  converged, 
solution  of  eigenvalues  and  the  corresponding  eigenfunction.  The  case  that  will  be 
examined  corresponds  to  6  =  60  degrees  in  Table  (4.3)  obtained  from  Carpenter 
and  Morris  [23)  with  R  =  2240,  and  u 7  =  0.055.  Carpenter  and  Morris  chose  an 
appropriate  compliant  coating  density  since  rubber  may  have  a  density  of  960-1300 
kg/m3  [41].  A  swivel-arm  angle  of  0  =  60  for  the  present  calculations  enables  a 
comparison  to  be  made  with  the  results  obtained  from  the  sixth-order  model  of 
Carpenter  and  Morris.  The  complex  wavenumber  indicating  the  onset  of  instability 


f  v  «- 


is  of  interest  in  this  investigation,  so  the  convergence  and  accuracy  of  such  are 
computed.  In  Table  (4.4)  the  wavenumber  is  shown  to  converge,  but  a  large  number 
of  polynomials  are  required  for  a  desired  accuracy.  The  rigid  wall  case  requires  only 
about  one-third  as  many  polynomials  for  a  comparable  accuracy.  Carpenter  and 
Morris  chose  48  polynomials  for  their  stability  calculations.  With  this  choice  the 
two-digits  of  accuracy  obtained  were  sufficient  for  obtaining  adequate  results.  The 
cost  of  additional  accuracy  may  be  seen  in  Table  (4.4)  where  a  gain  of  two  significant 
digits  results  in  approximately  triple  the  computational  cost. 

By  looking  at  the  least  damped  wavenumber  for  the  TSI  wave  over  a  frequency 
variation,  the  frequency  at  which  the  largest  growth  rate  occurs  may  be  determined. 
Shown  in  Figures  (4.6)  and  (4.7)  are  plots  of  the  wavenumber  verses  the  frequency 


Table  4.4:  Number  of  Chebyshev  polynomials  required  for  eigenvalue  convergence 
for  R=2240,  uJ  =  0.055,  0  =  60  and  B=0.08673  x  10~G. 


N 

a 

cpu  time(s) 

32 

.15805932, -.30840550  x  10_i 

40 

.15799165, -.31577908  x  10~2 

48 

. 15781832, -.31367360  x  10~2 

31.9 

52 

. 15780772,-. 31282229  x  10"2 

56 

.15780989, -.31412738  x  10-2 

60 

. 15781560,-. 31414286  x  10“2 

62 

. 15781491,-. 31393189  x  10~2 

65 

. 15781568,-. 31395635  x  10"2 

68 

.15781581, -.31399292  x  10~2 

69 

.15781518, -.31392974  x  10"2 

72 

.15781542, -.31396316  x  10“2 

75 

.15781538, -.31396114  x  10~2 

109.07 
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Figure  4.6:  Imaginary  part  of  the  complex  wavenumber  plotted  against  fre¬ 
quency  for  R=2240. 

—  TSI  from  [23]  O-present  calculations 
-  -  Traveling  wave-flutter  from  1 2 3 j  X-present  calculations 
1-9  =  0  2-9  =  60. 
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compared  with  the  results  obtained  by  Carpenter  and  Morris.  The  results  are 
in  agreement.  These  figures  serve  two  basic  functions.  First,  the  revised  model 
formulation  by  Morris  which  is  being  used  in  the  present  calculations  is  shown  to  be 
adequate  in  comparison  to  the  higher  order  model  of  [23).  And  second,  they,  along 
with  the  surface  properties  in  Table  (4.3)  form  a  starting  point,  or  reference  point,  in 
the  optimization  procedure.  In  conjunction  with  this,  a  comparison  between  these 
results  with  a  rigid  surface  makes  evident  the  possibility  of  delaying  transition.  At  a 
frequency  of  0.055  the  compliant  surface  has  a  least  damped  wavenumber  of  -.0031 
while  the  rigid  surface  has  a  value  of  -.01.  This  holds  for  the  surface  properties  of 
9  =  60.  But  this  investigation  seeks  to  show  that  by  varying  the  surface  properties, 
reduced  growth  rates  of  instability  or  even  complete  stabiliz?  +  ion  of  the  boundary 
layer  theoretically  may  be  achieved. 

A  question  as  to  why  the  large  number  of  Chebyshev  polynomials  is  required 
arises  for  the  compliant  case;  one  possible  answer  may  be  found  from  analyzing 
the  eigenfunction  behavior.  The  eigenfunctions  for  N=10,24,  and  48  are  shown  in 
Figure  (4.8)  for  the  least  damped  wavenumber  of  TSI.  The  corresponding  results 
show  rather  significant  differences  between  the  curves  near  the  boundary.  If  one 
were  to  make  a  comparison  of  the  numerical  aspects  between  the  rigid  and  compliant 
cases,  more  insight  may  be  shed  on  the  problem  in  question.  The  Chebyshev  series 
coefficients  in  general  have  the  property  that  the  leading  coefficient  is  the  largest 
in  magnitude.  The  remaining  coefficients  progressively  get  smaller  as  the  order  of 
the  terms  increase.  With  such  the  very  small,  normally  insignificant,  trailing  terms 
may  be  neglected  to  obtain  an  accurate  solution.  If  this  were  not  the  case,  then 
essentially  an  infinite  number  of  terms  would  be  required  for  a  solution.  For  the 
rigid  wall  case,  the  smaller  terms  may  be  neglected  and  a  relatively  accurate  solution 
is  obtained.  The  compliant  case  behaves  in  an  unconventional  manner.  The 


leading  four  or  five  coefficients  decrease  in  magnitude  gradually  as  is  expected;  the 
remaining  coefficients  drop-off  to  small  values  very  rapidly  and  non-uniformly.  As 
before,  it  might  be  expected  that  the  smaller  trailing  terms  may  be  dropped  and  an 
accurate  solution  would  be  achieved  requiring  fewer  Chebyshev  polynomials.  This  is 
a  somewhat  true  statement  since  for  N=10  a  rather  crude  approximation  is  achieved. 
On  the  other  hand,  a  sufficiently  accurate  solution  requires  the  very  small  trailing 
coefficients  to  remain  a  part  of  the  solution.  A  possible  reason  for  the  necessity  of 
the  additional  terms  may  lie  in  the  convergence  characteristics  of  the  function.  It 
may  be  possible  that  although  the  Euclidean  norm  of  the  system  becomes  small, 
or  converges,  this  may  not  be  a  sufficient  convergence  criteria.  Rather  the  infinity- 
norm  may  not  be  small. 

As  a  means  to  reduce  the  required  number  of  Chebyshev  polynomials,  stretch¬ 
ing  factors  were  implemented  to  decrease  the  amount  of  stiffness  in  the  problem. 
No  reduction  in  the  number  of  required  Chebyshev  polynomials  resulted.  Instead  of 
using  Chebyshev  polynomials  for  the  series  approximation,  improved  convergence 
of  the  series  might  be  realized  by  using  a  different  polynomial  such  as  the  Legendre 
polynomial.  Alternatively,  a  multi-domain  approach  may  be  attempted.  The  first 
domain  would  extend  from  the  compliant  surface  out  in  the  cross-stream  direction 
a  small  distance.  The  outer  domain  is  matched  with  the  inner  and  proceeds  to 
infinity.  The  solution  of  the  inner  domain  would  require  a  larger  number  of  poly¬ 
nomials  for  an  accurate  solution  as  compared  to  the  outer.  The  idea  behind  such 
an  approach  is  that  the  combined  solution  may  require  less  polynomials  than  the 
single  domain  problem. 

A  complete  explanation  for  the  behavior  of  the  series  approximation  for  the 
compliant  surface  problem  is  unconclusive  at  present. 
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CHAPTER  5 

EIGENMODE  SENSITIVITY  TO  BOUNDARY  PARAMETERS 

5.1  Introduction 

An  aspect  of  the  surface  property  optimization  is  obtaining  a  technique  to 
determine  the  changes  in  the  eigenmodes  with  respect  to  the  boundary  parameters. 
A  method  which  appears  to  have  potential  was  used  by  Bridges  and  Morris  [27,33] 
to  determine  the  frequency  of  the  most  unstable  eigenvalue  will  be  used  here  to 
determine  the  sensitivity  of  the  least  damped  eigenmodes  to  boundary  parameter 
changes.  In  Chapter  3  it  was  mentioned  that  Lancaster’s  local  eigenvalue  refinement 
method  could  be  extended  to  perform  a  portion  of  the  optimization  procedure.  The 
formulation  of  the  method  will  be  described  and  tested  using  the  model  problem 
then  extended  to  the  physical  problem. 

5.2  Model  Problem  Parameter 

Since  the  optimization  desired  occurrs  with  respect  to  boundary  parameters, 
a  modification  is  made  to  the  model  problem  (3.32).  A  nondimensional  surface 
parameter,  d,  is  introduced  giving  the  modified  boundary  condition. 

0aU(- 1)  +0'(-l)  =0  (5.1) 

The  spectral  discretization  results  in 

[D3(a^)|{a}  =  {0},  (5.2) 

where  {a}  is  the  right  eigenvector  and  the  lambda  matrix  is  given  by 


lal.r 


4, 
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D3{a,0)  =  Co0a3  +  Cja2  +  C2a  +  C3. 

Before  deriving  the  necessary  relation,  an  additional  vector  must  be  defined  and 
a  means  of  obtaining  this  vector  as  well  as  the  right  eigenvector  must  be  discussed. 
This  vector  may  be  defined  such  that 

{<T}H[D3(a,/?)]  =  {0},  (5.3) 

where  {a"}  is  referred  to  as  the  left  eigenvector  and  H  denotes  the  complex  conju¬ 
gate  transpose.  Relation  (5.3)  may  also  be  cast  in  the  form 

|D3"(5,(J)|{a-}  =  {0}.  (5.4) 

This  has  a  similar  form  to  (5.2)  for  the  right  eigenvector;so  a  common  technique 
for  determining  the  eigenvectors  may  be  used.  To  compute  a  single  eigenvector  the 
following  inverse  iteration  is  used 

|C(5)|{i‘+1}  =t{x‘},  (5.5) 

where  o  is  a  scaling,  or  normalizing,  factor.  The  procedure  converges  in  two  or  three 
iterations  using  an  initial  guess  of  {x0}  =  [1,1,  ...,l]T.  The  right  eigenvector  may  be 
conveniently  computed  using  the  already  formed  lambda  matrix,  [£>),  from  the  local 


method;the  left  eigenvector  may  similiarly  be  computed  with  the  inverted  lambda 
matrix,  [Z?_1j,  from  the  local  method  by-way-of  the  relation  [/1_1]/7  =  [/4W]~'. 

The  necessary  components  for  the  differentiation  have  been  computed,  so  the 
derivation  of  the  sensitivity  relations  will  follow.  If  the  matrix  system  (5.2)  is 
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,  1*  *  • 


differentiated  with  respect  to  the  boundary  parameter,  (3,  the  following  may  be 
obtained. 


+  +  ={0}  (5.6) 

By  multiplying  (5.6)  by  {a'}H  the  last  term  is  eliminated  and  the  result  is 

^  _  _  {a'}  (d[D3\/d/3)  {a} 

3D  {a-)»[0<"(a,fl)|{a)' 

Recall  that  the  parameter,  /?,  appears  only  in  the  leading  coefficient  matrix,  so  (5.7) 
may  be  given  by 


^  = _ {a'}tf;C0]a3{a} _ 

df3  {q‘}"[3/?C0Q2  +  2C,a  +  C2){q}'  [  ’ 

The  matrix  [Z?^]  may  a*so  be  taken  from  the  local  method  described  in  Chapter 
3.  From  (5.8)  a  means  has  been  obtained  to  determine  the  effect  of  an  eigenmode 
to  changing  surface  parameters,  or  more  specifically,  the  sensitivity  of  the  least 
damped  wavenumber  to  surface  property  variations. 

The  accuracy  of  (5.8)  may  be  determined  by  a  comparison  of  this  method 
with  a  finite  difference  approximation.  The  results  in  Table  (5.1)  for  N  =  11  and 
Zj  =  0.25  show  good  agreement  between  the  finite  difference  approximation  and  the 
approach  of  (5.8)  for  the  model  problem. 

5.3  Compliant  Surface  Parameters 

The  formulation  for  the  sensitivity  of  the  eigenmode  to  a  boundary  parameter 
may  similiarly  be  applied  to  the  mechanical  model  representing  a  compliant  surface. 
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Table  5.1:  Sensitivity  of  the  eigenvalue,  a,  to  the  surface  parameter,  /3,  in  the 
model  problem  with  N=ll  and  cJ  =  0.25. 


R 

0 

a 

Numerical 

F.Diff. 

Error 

1. 

1.00 

.5454 

-.8344  x  10-1 

-.8492  x  10_1 

.17  x  10"3 

0.95 

.5496 

-.8642  x  10~ 1 

-.8645  x  10" 1 

.33  x  10-° 

0.90 

.5540 

-.8960  x  10_1 

-.8799  x  10" 1 

.18  x  10"3 

100. 

1.00 

.9451 

-.2756  x  10_1 

-.2779  x  10_1 

.12  x  10“3 

0.95 

.9465 

-.2821  x  10-1 

-.2822  x  10_1 

.16  x  10"5 

0.90 

.9479 

-.2890  x  10" 1 

-.2855  x  10-1 

.12  x  lO-3 

500. 

1.00 

.09953 

-.2194  x  10~3 

-.2195  x  10-3 

.35  x  10"5 

0.95 

.09954 

-.2196  x  10-3 

-.2196  x  10~3 

.34  x  10“° 

0.90 

.09955 

-.2197  x  10"3 

-.2197  x  10~3 

.29  x  10"5 

v\> 


■  V  V. 
C'-VA 
■*>.> 
/  ".*  ■> 

>■  V  •  “  • 


v 

.v.v.v 

V  V  Vi 


(5.116) 


(«'}"  [(o=36^trt-g)/(p<,C/^3f  ){v'(  — 1)}|  {a} 

{<■■}«  [£>S,,(5)I  fo> 

{<■'}"  [(5<-)/(^t/^,){g'(-l)}|  {a} 

{a'}»  |Cj1)(5)|  {a) 

K}"  [(-^2l>)/(g^'){t',(-l)}l  W 
{' «■}"  [^‘>(5)1  M 


(5.11c) 


(5. lid) 


(5. lie) 


Both  approaches  should  render  identical  results.  The  reason  for  giving  both  ap¬ 
proaches  is  to  verify  the  prospect  of  performing  the  analysis  with  the  four  nondi- 
mensional  variables  as  opposed  to  the  five  physical  variables.  Depending  on  the 
scope  of  the  analysis  for  a  given  problem  using  the  above  differential  technique,  it 
may  be  more  practical  to  use  fewer  nondimensional  variables  if  possible. 

As  with  the  model,  the  method  is  shown  to  be  accurate  in  comparison  to 
finite  difference  approximations.  The  results  are  compared  in  Table  (5.2)  with 
N  =  48,  u7  —  0.055,  R  =  2240,  and  the  surface  properties  corresponding  to  the  case 
where  0  =  60  of  Carpenter  and  Morris  [23].  A  comparison  showing  the  influence  of 
the  number  of  Chebyshev  polynomials  on  the  accuracy  of  the  differential  is  given  in 
Table  (5.3).  Obviously  increasing  the  number  of  Chebyshev  series  terms  increases 
the  number  of  accurate  significant  digits.  To  obtain  an  eigenvalue  with  two  sig¬ 
nificant  digits  of  accuracy  48  polynomials  are  required.  For  the  sensitivity  with  a 
similiar  number  of  polynomials  two  significant  digits  of  accuracy  are  retained.  This 
may  lead  one  to  conclude  that  the  accuracy  of  the  sensitivity  measure  is  related 
to  the  eigenvalue  accuracy  of  interest.  So  the  choice  of  the  number  of  Chebyshev 
polynomials  required  is  based  on  the  desired  accuracy  of  the  eigenvalue  in  question. 

Up  to  this  point  the  tools  necessary  for  the  optimization  procedure  have  been 
derived  and  analyzed.  The  local  eigenvalue  refinement  method,  the  eigenvector 
determination  scheme,  and  the  form*  la  relating  the  sensitivity  of  eigenmodes  to 
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Table  5.3:  Variation  of  the  sensitivity  of  the  eigenmode  to  changes  in  the  boundary 
property,  B,  with  the  number  of  Chebyshev  polynomials.  (B{Nm)  =  0.773  x  10-7). 


W 

M 

m 


.15804017, -.00293839 
.15757249, -.00295601 
.15774825, -.00295831 
.15781629, -.00313296 
.15781285, -.00314088 


Numer:da/3B 

213.16244,-78.565489 

217.97623,-73.151999 

216.67877,-75.295445 

215.11468,-74.643011 

214.94787,-74.624582 


F.DifT.:da/dfl 

101.4,-80.12 

218.0,-74.39 

216.0,-76.20 

215.0,-74.72 

215.0,-74.64 


CHAPTER  6 


OPTIMIZATION  OF  BOUNDARY  PARAMETERS 

6.1  Introduction 

In  the  remaining  segment  of  this  discussion  the  importance  of  the  surface  prop¬ 
erty  changes  on  boundary  layer  stabilization  will  be  shown.  Numerical  techniques 
outlined  in  the  previous  chapter  are  utilized  to  measure  the  instability  growth  or 
decay  depending  on  the  surface  property  variations. 

6.2  Minimization  Method  and  Results 

In  Chapter  5  a  means  was  developed  to  determine  numerically  the  sensitivity 
of  an  eigenvalue  to  surface  property  changes.  Relatively  speaking,  this  may  be 
considered  a  gradient.  This  is  verified  by  a  comparison  with  the  finite  difference 
approximation.  These  gradient  measurements  are  the  desired  agents  needed  for  the 
optimization  algorithm  available. 

Many  algorithms  have  been  proposed  for  minimizing  a  function  of  multi- 
variable  problems.  Gradients  are  not  necessary  to  obtain  the  desired  results  as 
is  explained  by  Powell  [42j;  however,  many  benefits  are  found  in  the  use  of  gradi¬ 
ents.  For  instance,  the  relative  influence  of  each  property  is  made  evident  in  the 
gradient.  In  the  present  investigation,  the  influence  of  a  property  may  be  found  to 
have  a  dominant  effect  on  the  growth  or  decay  of  the  instability  as  opposed  to  a 
property  having  no  effect  on  the  instability.  This  is  significant  in  the  present  prob¬ 
lem  for  understanding  the  relative  importance  of  the  sensitivity  of  TSI  and  FISI 
waves  with  the  surface  property  changes.  The  measure  of  TSI  sensitivity  was  given 
in  Table  (5.2),  and  the  FISI  sensitivity  measurements  are  given  in  Table  (6.1).  The 
TSI  sensitivity  values  are  one  or  two  orders  of  magnitude  greater  than  the  FISI 


Table  6.1:  Sensitivity  of  the  imaginary  part  of  the  wavenumber  of  FISI  to 
compliant  surface  property  changes  for  6  =  60  (Carpenter  and  Morris)  and 
at  =  .1462  x  10"3. 


da/ db(mm) 

da/dpm(kg/m3) 

da/dB(Nmm) 

da/dE(N/mm2) 

da/dK{N/mm3) 


+  .14149130  x  10~3 
+  .23184345  x  10“7 
-.24979009  x  10"5 
-.14693146  x  10~4 
-.21739922  x  10"2 
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values.  In  theory,  the  stabilization  of  one  class  of  instability  leads  to  a  destabiliza¬ 
tion  of  the  other  class.  From  the  results  at  present,  a  means  has  been  developed 
to  measure  the  relative  stabilizing-destabilizing  effect  of  the  two  classes  of  instabil¬ 
ities.  This  is  significant  since  it  indicates  the  dominant  influence  of  TSI  waves  in 
reference  to  the  choice  of  compliant  surface  properties.  So  the  optimization  should 
be  concerned  with  the  stabilization  of  the  more  critical  TSI  waves. 

From  this  comparison,  a  multi-variable  algorithm  is  sought  to  minimize  the 
TSI  growth  rate.  Such  an  algorithm  was  made  available  by  Dr.  William  Hagar  of 
the  Mathematics  Department  at  Penn  State.  The  algorithm  is  based  on  the  con¬ 
jugate  gradient  method  by  Fletcher  and  Reeves  [43 j .  The  code  routinely  calls  the 
surface  properties,  gradients,  and  least  damped  wavenumber  during  the  iterations. 
Beginning  with  the  properties  of  Carpenter  and  Morris  for  9  =  60,  iterations  are 
performed  and  listed  in  Table  (6.2).  For  three  iterations  a  stabilization  begins  to 
occur.  This  trend  would  continue  and  lead  to  a  damped  wavenumber.  For  each  it¬ 
eration  a  systematic  repetition  of  property  values  is  observed.  This  may  occur  due 
to  the  relati  ve  invarience  of  the  sensitivity  values.  In  using  this  algorithm  the  prop¬ 
erties  must  have  comparable  magnitudes  or  the  iterative  process  leads  to  physically 
unrealizable  values.  For  example,  the  flexural  rigidity  which  has  a  small  magnitude 
was  varied  and  became  negative.  This  may  be  avoided  by  limiting  the  band  of 
possible  property  values  available  to  the  algorithm.  This  is  unnecessary  at  present 
since  it  proves  more  efficient  computationally  to  use  a  simple  variational  approach. 
Although  exact  cpu  accounts  are  unavailable,  a  relative  comparison  is  possible.  The 
above  method  requires  appoximately  five  minutes  on  the  VAX  11/8580  as  opposed 
to  four  minutes  for  the  approach  that  will  follow.  The  above  method  did  confirm 
that  stabilization  is  theoretically  possible  through  the  appropriate  surface  property 
combinations,  but  much  more  useful  information  is  obtained  from  the 


Table  6.2:  Minimization  of  instability  growth  rate  by  the  conjugate  gradient  ap¬ 
proach  for  B(Nm)=0.773  x  10-7.  pm{kg j m3)  =  1000,  and  an  initial  step  of  0.05. 


Iteration 

b 

(mm) 

E 

(N/mm2) 

K 

(N/mm3) 

a, 

Initial 

.11100 

.50900 

.05900 

-.003133 

I 

.10836 

.50835 

.05559 

-.002791 

.10953 

.50864 

.05740 

-.002946 

.10836 

.50835 

.05559 

-.002791 

2 

.10716 

.50806 

.05402 

-.002622 

.10769 

.50819 

.05472 

-.002698 

.10716 

.50806 

.05402 

-.002622 

3 

.10679 

.50797 

.05351 

-.002567 

.10695 

.50801 

.05373 

-.002592 

.10679 

.50797 

.05351 

-.002567 

method  that  will  follow. 


6.3  Variational  Method  and  Results 


For  multi-variable  problems  it  is  advantageous  to  fix  some  variables.  The  char¬ 
acteristics  of  the  remaining  variables  in  the  problem  may  then  be  observed.  In  the 
present  discussion  the  flexural  rigidity,  thickness,  and  modulus  of  elasticity  of  the 
plate  are  functionally  related,  so  the  remaining  parameters  are  fixed.  This  naturally 
occurs  for  the  plate  density  since  the  sensitivity  measurements  in  Table  (5.2)  show 
that  the  instability  is  not  influenced  significantly  by  density  changes. 

In  this  analysis  a  simple  variation  of  properties  is  made.  The  properties  are 
governed  by  the  following  relation 


12(1  -  t/2)’ 


where  the  Poisson’s  ratio,  u,  is  0.5.  Two  approaches  were  used.  The  first  approach 
maintains  an  essentially  constant  value  of  the  flexural  rigidity  and  appropriately 
varies  the  thickness  and  modulus  of  elasticity  of  the  plate.  The  results  are  given  in 
Table  (6.3)  and  plotted  in  Figures  (6.1)  and  (6.2)  against  the  least  damped  wave- 
number  of  TSI.  The  results  indicate  that  by  simultaneously  increasing  the  plate 
thickness  and  decreasing  the  modulus  of  elasticity  a  stabilization  of  the  boundary 
layer  may  be  realized.  The  second  approach  maintains  essentially  constant  values 
of  the  thickness  and  modulus  of  elasticity  and  varies  the  flexural  rigidity.  These 
results  are  given  in  Table  (6.4)  and  plotted  in  Figures  (6.3)  and  (6.4).  These  results 
indicate  that  by  decreasing  the  flexural  rigidity  the  boundary  layer  tends  to  become 
stabilized.  The  most  pronounced  stabilization  occurs  when  the  plate  thickness  is 
increased  and  the  modulus  of  elasticity  is  decreased.  In  both  cases  the  FISI  values 
show  little  change.  This  may  be  expected  as  indicated  by  the  sensitivity 
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.1985 
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.1458x10 

.1472x10 

.1486x10 

.1500x10 

.1512x10 

.1526x10 

.1540x10 

.1553x10 

.1567x10 

.1579x10 

.1592x10 

.1606x10 

.1620x10 

.1633x10 

.1647x10 
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.1675x10 
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instability  versus  15  and  b  where  b  is  primarily  chanced 
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Table  6.4:  Sensitivity  of  least  stable  wave-number  of  TSI  and  FISI  to  changes  the 
surface  properties:  B,E,  and  b  with  K(N/mm3)  =  .059  and  pm(kg/m3)  =  1000. 


(Nm) 

(mm) 

(N/mm2) 

a, 

a. 

.0743x10"° 

.1035 

.6031 

—  .3731xl0-2 

,1443xl0-3 

.0763x10"° 

.1035 

.6194 

—  .3878xl0-2 

.1441xl0"3 

.0783x10"° 

.1035 

.6356 

-.4020xl0-2 

.  1439xl0"3 

.0803x10-° 

.1035 

.6518 

-,4158xl0-2 

,1437xl0"3 

.0823x10-° 

.1035 

.6681 

-.4289xl0-2 

.1435xl0"3 

.0843x10-° 

.1085 

.5940 

—  .3917xl0-2 

.1448xl0"3 

.0863x10-° 

.1085 

.6081 

-.4014xl0-2 

.1446xl0"3 

.0883x10"° 

.1085 

.6222 

—  .4139xl0"2 

.1445xl0~3 

.0903x10-° 

.1085 

.6363 

-,4262xl0-2 

.1443xl0_3 

.0923x10-° 

.1085 

.6504 

-.4381xl0-2 

.1441xl0"3 

.0943x10"° 

.1135 

.5805 

—  .4016xl0"2 

.1455xl0"3 

.0963x10'° 

.1135 

.5928 

-.4102xl0-2 

.1453xl0"3 

.0983x10-° 

.1135 

.6051 

-.4214xl0-2 

.  1452xl0"3 

.1003x10-° 

.1135 

.6174 

. 

—  .4324xl0"2 

.  1450xl0"3 

.1023x10-° 

.1135 

.6297 

-.4431xl0-2 

. 1448xl0-3 

.1043x10-° 

.1185 

.5641 

—  .4074xl0"2 

.1462xl0"3 

.1063x10"° 

.1185 

.5749 

-  .4 15  IxlO"2 

.1461xl0"3 

.1083x10-° 

.1185 

.5858 

-  4253cl0"2 

.1459xl0-3 

.1103x10-° 

.1185 

.5966 

—  .4353xl0"2 

,1457xl0"3 

.1123x10"° 

.1185 

.6074 

—  .4451xl0"2 

.1456xl0"3 

.1143x10-° 

.1235 

.5461 

—  .4101xl0"2 

.1470xl0"3 

.1163x10"° 

.1235 

.5557 

j 

-.4171xl0-2 

.  1468xl0"3 

.1183x10"° 

.1235 

.5652 

-.4265xl0_2 

.1467xl0"3 

.1203x10-° 

.1235 

.5748 

-  .4357xl0-2 

.1465xl0"3 

.1223x10"° 

.1235 

.5843 

-  .4448xl0-2 

.1464xl0-3 

-0.44d-2 
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measurements.  A  more  indepth  explanation  may  be  desired  to  justify  that  the  re¬ 
sults  obtained  with  respect  to  FISI  growth  rate  invariation  are  in  agreement  with 
theory.  Carpenter  and  Garrad  [18]  developed  a  means  to  identify  curves  indicat¬ 
ing  expected  stable  and  unstable  regions  for  FISI.  These  curves  are  dependent  on 
the  properties  of  an  isotropic  surface  for  the  temporal  case.  If  their  analysis  were 
extended  to  the  spatial  non-isotropic  case,  it  seems  reasonable  to  expect  that  “sta¬ 
bility  curves”  may  be  determined.  The  property  variation  in  the  present  analysis 
may  ensure  that  FISI  remains  in  a  neutral  or  stable  region.  It  may  be  found  that 
the  destabilization  of  FISI  waves  occurring  due  to  an  increase  in  the  plate  thickness 
may  be  offset  by  the  corresponding  stabilization  occurring  due  to  an  increase  in  the 
flexural  rigidity. 

The  final  variation  is  made  with  respect  to  the  model  swivel-arm  angle.  Shown 
in  Figure  (6.5)  is  the  effect  on  the  least  damped  wavenumber  for  TSI  to  changes  in 
the  swivel-arm  angle  for  fixed  surface  properties.  These  properties  correspond  to 
those  given  by  Carpenter  and  Morris  for  6  -  60  with  R  —  2240  and  uJ  =  0.055.  As 
was  found  with  the  surface  property  changes,  a  stabilization  may  be  realized.  This 
can  primarily  be  achieved  with  an  angle  between  0  and  50.  It  should  be  noted  that 
the  isotropic  case  corresponding  to  6  -  0  has  a  stabilizing  affect  on  the  boundary 
layer  for  these  particular  surface  properties.  Reynolds  number  and  frequency.  This 
concludes  t  he  findings  discovered  through  this  investigation.  Of  course,  further  cases 
may  be  performed,  but  this  example  shows  the  relative  importance  of  the  surface 
properties  on  the  instability  and  properties  which  may  lead  to  a  stabilization  of  the 
boundary  laver. 


0.56d-2 


Variation  of  the  least  stable  wave-number  of  Tollmien-Schlichting 


CHAPTER  7 


DISCUSSION  AND  CONCLUSIONS 

The  majority  of  this  thesis  has  consisted  of  a  description  of  the  problem  and 
the  building  of  the  numerical  tools  necessary  to  fulfill  the  expectations  of  the  inves¬ 
tigation.  The  technique  to  determine  the  sensitivity  of  the  instability  with  respect 
to  surface  property  changes  was  critical  in  understanding  the  importance  of  each 
property.  This  measurement  also  provided  a  means  to  measure  the  relative  influence 
of  the  surface  properties  on  both  TSI  and  FISI  waves.  The  minimization  algorithm 
seeks  a  decreasing  gradient  as  a  convergence  criterion.  In  this  analysis,  however,  the 
gradients  changed  very  little  over  the  property  range  of  interest.  So  an  algorithm 
which  travels  in  a  gradient  dependent  manner  is  less  practical  than  a  simple  prop¬ 
erty  variation  for  the  present  problem  of  interest.  It  also  proves  computationally 
more  efficient  and  much  more  information  is  gained  when  the  property  variation 
approach  is  used.  In  extending  the  present  analysis,  the  surface  properties  may 
be  sought  which  give  minimal  sensitivity  values.  If  this  were  attained,  one  might 
expect  that  with  small  changes  in  flow  conditions  the  instability  growth  rate  would 
essentially  remain  unchanged.  The  idea  seems  plausible  but  in  fact  the  measure  of 
sensitivity  over  the  range  of  present  property  values  did  not  change  significantly  in 
magnitude.  This  may  be  shown  graphically  by  reviewing  the  slopes  of  the  curves 
in  Figures  (6.1  )-(6.4) . 

In  the  final  method,  the  variation  approach  gives  a  simple  means  to  attain 
stability  curves  with  respect  to  surface  properties.  In  a  more  complete  sense,  for  a 
multi-variable  problem  it  is  possible  to  attain  “stability  planes-’ .  In  this  investiga¬ 
tion  the  range  of  surface  properties  centered  around  the  values  by  Carpenter  and 
Morris.  If  manufacturable  surface  property  combinations  were  available,  it  would 


be  possible  to  predict  a  coating  most  likely  to  delay  transition.  Assuming  Kramer’s 
conclusions  concerning  surface  imperfections  and  water  impurities  resulting  in  no 
performance  loss  hold,  the  predicted  performance  should  be  realized  if  the  manu¬ 
factured  compliant  coating  is  in  accord  with  the  mechanical  model  representation. 

From  the  results  obtained  in  this  investigation,  further  research  may  commence 
in  many  directions.  Experimentally,  a  surface  may  be  constructed  and  tested  on 
a  model.  The  performance  may  then  be  compared  with  an  uncoated  model  and 
the  predicted  results  of  the  mechanical  model.  Along  a  similiar  route  taken  in  this 
study,  the  delay  to  transition  may  be  analyzed  for  an  optimal  set  of  surface  prop¬ 
erties  over  a  range  of  flow  conditions.  This  is  an  obvious  necessity  for  commercial 
considerations.  Of  course,  compliant  coatings  are  not  limited  to  laminar  transition 
analysis  and  are  also  being  used  in  turbulence  research.  It  would  be  of  interest  to 
determine  the  desirable  compliant  coating  properties  in  turbulent  flow.  These  prop¬ 
erties  could  be  compared  with  the  “stability  planes”  which  may  be  obtained  from 
the  present  analysis.  Overlap  regions  may  be  found  which  when  used  in  practice 
delay  transition  and  perform  favorably  after  the  transition  to  turbulent  flow.  This 
would  enhance  the  performance  of  a  coated  vehicle  over  a  range  of  velocities. 


APPENDIX  A 


CHEBYSHEV  SERIES  FORMULAE 

The  Chebyshev  polynomials,  Tn(x),  are  defined  on  the  interval  x  6  (-1,  +  l| 
and  are  derived  from  and  related  to  the  cosine  function  by 


Tn  (cosff)  —  cosnd,  (A.l) 

with  the  initial  few  polynomials  appearing  as  Ta( x)  =  l,Ti(x)  =  x, 

T2 (x)  =  2x2  -  l,r3(x)  =  4x3  -  3x,  ect.  The  following  trigonometric  identity  can  be 
obtained. 


cos(n  +  1)0  =  2 cosO  ■  cosnO  -  cos(n  -  1)0 


(A.2) 


This  results  in  a  Chebyshev  recurrence  formula  for  higher  order  polynomials. 


rn  +  i(x)  =  2 xTn{x)  -  r„_i(x) 


The  product  formula  is 


1 


Tn{x)Tm{x)  —  ^  (Tn  +  m(x)  +  T|n_m| (x)) 


and  the  indefinite  integral  relation  is 


T  I 


r)  r/  -r  — 


Tt(x) 


(A.3) 


(A.  4) 


1  (To(x)  +  T2(x)) 


n  =  0 
n  —  1 


(  J 


The  series  boundary  conditions  are 


r„(±l)  =  (±l)n  (.4.6) 

and  the  differential  relation  for  Chebyshev  polynomials  at  the  boundaries  is 

jp  p~1 

— r„(±i)  =  (±  rtf"2  -  k7VVk  +  »•  M-7) 

k  =  0 

Another  efficient  relation  useful  when  performing  the  summation  of  a  Cheby¬ 
shev  series  is  given  by 


Y  ’anTn{x)  =  l~  [6„(i)  -  62(x)j  ,  (.4.8) 

n~o 

where  the  prime  signifies  that  the  leading  term  is  to  be  halved.  The  recurrence 
system  needed  to  evaluate  (A. 8)  is 

bn{x)  =  2x6n  +  ,(x)  -  bn  +  2(x)  +  an  (A. 9) 

6/v  +  i(x)  =  by  4. 2  ( x )  =  0. 

A  Chebyshev  formula  useful  in  approximating  a  known  function  in  a  Chebyshev 
series  can  be  defined  as 


,v 

*(*)  =  Y'^Tn(x),  (.4.10a) 

n  =0 


where  <t>(x)  is  a  known  function  at  all  points  in  the  domain.  The  coefficients.  on  . 
are  given  by 


(A. 106) 


with 


fc=0 


Xk 


cos 


krr 

iV 


for  k  =  0,1,2, 


(A. 10c) 


The  double  prime  on  the  summation  signifies  that  the  leading  and  trailing  coeffi¬ 
cients  are  to  be  halved.  The  final  Chebyshev  property  that  will  be  given  prior  to 
listing  practical  integral  formulae  is  the  approximation  of  the  differential  of  a  known 
function  in  Chebyshev  series.  The  derivative  is  given  by 


OO 

0'(x)  =  ^  6nr„(x)  (A. 11a) 

n  =  0 

where 


(A. 116) 


and 

(  2  n  =  0 

Cn  “  l  1  n  >  0. 

The  coefficients,  an  ,  are  obtained  from  the  series  approximation  to  the  known 
function, 

To  obtain  the  solution  of  a  differential  equation  by  a  Chebyshev  series  approxi¬ 
mation,  it  is  convenient,  although  not  necessary  to  convert  the  differential  equation 
to  an  integral  form.  As  such,  a  function  is  represented  by  the  following  finite, 
Chebyshev  series. 


(A. 11c) 


6„  = 


E 

p  =  n+  1 
p  +  nodd 


paP 
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4>{x)  =  ^  'anTn[x) 


(.4.12) 


By  applying  the  integral  relation  (A. 5)  appropriately  and  repeatedly,  the  following 


relations  are  obtained. 


where 


where 


.  N  +  i 

1.  /  (f>{x)dx  =  '  bnTni X) 

*  _ _ r\ 


bn  =  —  K-i  -  an+i)  for  n>  1 
In 


f  N  +  2 

2.  J  J  4>{x)dx 2  =  £  '  6nrn(x) 


L  _  an-2  an  ,  an  +  2  r  ^  „ 

6"  ~  [Z+TT)  -  i(^TT)  +  J^TT)  for  n-2 


f  f  r  'v+3 

3 ■  J  j  j  <t>(z)di3  =  Y2 '  bnTnW 


where 


(A. 13a) 


(A. 136) 


(A. 14a) 


(A. 146) 


(A. 15a) 


^  _  _ fln — 3 _  3an_ i  3an  + 1 

8n(n  -  l)(n  -  2)  8n(n  —  2)(n  +  1)  8n(n  —  l)(n  -t-  2) 

-  Z~7 - 777 - 7T  f°r  n  >  3  (A. 156) 

8n(n  +  1  n  +  2  _  v  ’ 


where 
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.  Q-n  —  4  On  —  2  n 

n  16n(n  —  l)(n  —  2 )(n  -  3)  4n(n2  -  l)(n  -  3)  8 (n2  —  l)(n2  -  4) 


4n(n2  —  l)(n  +  3)  16n(n  +  1  )(n  +  2  )(n  +  3) 


for  n  >  4  (A. 166) 


VVhen  the  coefficients  in  the  differential  equations  are  non-constant,  the  Cheby- 
shev  product  formula  (A. 4)  is  needed.  Introducing  a  function,  u(z),  representing 
the  non-constant  coefficient,  the  following  is  obtained. 


u{x)4>{x)  =  ^  'dnTn(x) 


(A. 17a) 


u(x)  =  Y2  'UnTn(x) 


(A. 176) 


1  1  v- 

dn  =  ~una0  +  -  2_^  (ulm-n[  +  Um+n)am  for  n  >  0 


(A. 17c) 


Integrations  are  performed  in  a  straight  forward  manner  using  the  integral  re¬ 
lation  (A. 5).  The  following  integral  relations  prove  useful  for  the  problem  presented 
in  this  thesis. 


r 

1.  /  u(x)<t>(x)dx  =  T  '  dnTn(x) 

*'  —  _ n 


(A. 18a) 


where 


dn  ( U  rt  -  l  Un4.[ 

4n 
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>  (”>m- 

4n  i — ' 

*r\  —  1 


n  1  m  -  n  -  l 


^  rr\  n  l  )  &  m  for  n  _•  1  (.1.186) 
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APPENDIX  B 

BLASIUS  SOLUTION 

REPRESENTED  BY  A  CHEBYSHEV  SERIES 

In  this  thesis,  the  solution  to  the  Blasius  equation  is  necessary  and  is  repre¬ 
sented  by  a  Chebyshev  series.  In  order  to  attain  this  end,  an  accurate  means  of 
obtaining  a  numerical  solution  is  first  necessary.  In  the  similarity  variables  the 
governing  equation  is  an  ordinary  differential  equation  of  the  form 

f'"  4-  -  f f"  =  0  (B. \) 

where  (  '  )  =  and  the  appropriate  boundary  conditions  are 

/(0)  =  f'(0)  --  0  .  lim  /'(>?)-  1.  (B.2) 

r;  — c 

where  the  derivative  of  /(rj)  is  the  streamwise  velocity.  Using  a  shooting-type 
method  with  a  5th-6th  order,  variable-step  Runge-Kutta  solver) IMSL:I)\  ERK)  . 
the  initial  condition  satisfying  the  streamwise  velocity  limit  can  be  found.  This 
results  in 

/"( 0)  -  0.3320573362185815  .  (B. 3) 

The  discrete  points  desired  in  the  ( 'hebyshev  domain,  z  •-  i  —  1 .  -e  1  j.  are  trans¬ 
formed  to  the  Blasius  variable,  rj  -  O.oc).  via  the  algebraic  transformation 

n  =  L  ■  (1  +  y ) / ( 1  -  y)  {B.4) 


where  y  is  the  physical  coordinate  with  domain  [0, oo).  The  solution  obtained  at 
the  desired  points  are  then  transformed  back  to  the  computational  domain  using 
the  inverse  of  ( B .-4 ) .  Taking  the  solution  of  a  functionfi.e.,  the  Blasius  solution), 
F{z),  given  at  all  points  in  z  6  |  —  1,  +  1).  a  Chebyshev  expansion  of  such  is  sought. 

,v 

F(z)  =  V'/nrn(2)  (B.  6) 

n =0 

The  prime  signifies  that  the  leading  term  of  the  series  is  to  be  halved.  An  exact 
solution  is  obtained  for  N  —  oo.  For  a  series  expansion,  the  function  must  be 
evaluated  at  the  Chebyshev  points 

z,  =  cos(ni/N)  i  =  0,1,2 .  {B. 7) 

The  series  at  these  points  is 

N 

F{zi)  =  53'/"r"(*‘)  ■  ,RS| 

n  =  0 

Using  the  relationship  between  the  Chebyshev  polynomial  and  the  cosine  funrtmn 
a  curve-fitting  formula  can  be  obtained.  Thus,  the  coefficients  can  be 
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where  the  double  prime  signifies  that  the  leading  and  trailing  terms  of  the  series 
are  to  be  halved.  With  the  identity 

Tn{zt)  =  cos(nni/ N)  (fl.10) 

the  desired  form  of  the  curve-fitting  formula  is  obtained.  Making  a  substitution  of 
(B.10)  into  (B.9),  the  following  results 

2  N 

fn  =  Jj  ^2"f(z,)c°S(nni/N)  .  (B.ll) 

1=0 

By  making  use  of  (B.ll)  with  (B.6),  the  Chebyshev  series  representation  of  a  func¬ 
tion  can  be  computed  to  a  desired  accuracy  by  taking  N  to  be  large. 

The  solution  of  the  Blasius  equation  represented  by  a  Chebyshev  series  is  at¬ 
tained  with  this  curve-fitting  formula.  Sufficient  accuracy  was  attained  using  an 
approximation  with  an  order  of  N  =  99.  The  solution  given  in  the  computational 
domain  is  shown  in  Figure  (3.1)  in  Chapter  3. 
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